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Single particle tracking in systems showing anomalous diffusion: the role of weak 

ergodicity breaking 
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Anomalous diffusion has been widely observed by single particle tracking microscopy in complex 
systems such as biological cells. The resulting time series are usually evaluated in terms of time av- 
erages. Often anomalous diffusion is connected with non-ergodic behaviour. In such cases the time 
averages remain random variables and hence irreproducible. Here we present a detailed analysis of 
the time averaged mean squared displacement for systems governed by anomalous diffusion, con- 
sidering both unconfined and restricted (corralled) motion. We discuss the behaviour of the time 
averaged mean squared displacement for two prominent stochastic processes, namely, continuous 
time random walks and fractional Brownian motion. We also study the distribution of the time 
averaged mean squared displacement around its ensemble mean, and show that this distribution 
preserves typical process characteristic even for short time series. Recently, velocity correlation 
functions were suggested to distinguish between these processes. We here present analytucal ex- 
pressions for the velocity correlation functions. Knowledge of the results presented here are expected 
to be relevant for the correct interpretation of single particle trajectory data in complex systems. 

PACS numbers: 02.50.-r, 05.40.Fb, 05.10.Gg 



I. INTRODUCTION 



Single particle tracking microscopy provides the pos- 
ition time series r(t) of individual particle trajectories 
in a medium [lHj]. The information garnered by single 
particle tracking yields insights into the mechanisms and 
forces, that drive or constrain the motion of the particle. 
An early example of systematic single particle trackingis 
given by the work of Jean Perrin on diffusive motion [5| . 
Due to the relatively short individual trajectories, Perrin 
used an ensemble average over many trajectories to ob- 
tain meaningful statistics. A few years later, Nordlund 
conceived a method to record much longer time series @ , 
allowing him to evaluate individual trajectories in terms 
of the time average and thus to avoid averages over not 
perfectly identical particles. Today, single particle track- 
ing has become a standard tool to characterise the micro- 
scopic rheological properties of a medium [7||, or to probe 
the active motion of biomolecular motors [8j. Particu- 
larly in biological cells and complex fluids single particle 
trajectory methods have become instrumental in uncov- 
ering deviations from normal Brownian motion of pass- 
ively moving particles 

Classical diffusion patterns are sketched in the left 
panel of Fig. [1] Accordingly, one may observe free diffu- 
sion, leading to a linear growth with time of the second 
moment [Line 1 in Fig. [Tj. Brownian motion may also 
be restricted (corralled, confined). Confinement in a cell, 
for instance, could be due to the cell walls. In that case 
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Figure 1: Diffusion modes of the time averaged mean squared 
displacement ((2]) as function of the lag time A. Left: Nor- 
mal diffusion growing like S 2 ~ A (1), restricted (confined) 
diffusion with a turnover from ~ A to ~ A (2), drift dif- 
fusion ~ A 2 (3). Right: Ergodic subdiffusion ~ A a (1'), 
restricted ergodic subdiffusion turning over from ~ A° to 
~ A (2'), non-ergodic subdiffusion ~ A (1"), restricted non- 
ergodic subdiffusion turning over from ~ A to ~ A 1_Q (2"), 
superdiffusion A 1+a (3'). Here, < a < 1. Note the double- 
logarithmic scale. 



the second moment initially grows linearly with time and 
eventually saturates to its thermal value equalling the 
second moment of the corresponding Boltzmann distri- 
bution [Line 2] . In the presence of a drift the second mo- 
ment grows with the square of time [Line 3] . Such results 
are typical for simple fluids. In more complex environ- 
ments, different patterns may be observed, as displayed 
on the right of Fig. [T] Here, subdiffusion may occur, 
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for which the second moment grows slower than linearly 
with time [Line 1']. Restricted subdiffusion would depart 
from this behaviour to reach a plateau [Line 2']. Driven 
motion may lead to a superdiffusive power-law form of 
the second moment with an exponent between 1 and 2 
[Line 3']. However, as we demonstrate below, subdiffu- 
sion may also be non-ergodic, and the associated time 
averaged second moment may grow linearly with time 
[Line 1"]. Similarly strange behaviour may be observed 
for restricted non-ergodic subdiffusion, which exhibits a 
power-law growth, not a saturation to a plateau [Line 
2" ] . Non-ergodic processes come along with a significant 
scatter between individual trajectories. This is an effect 
of the ageing nature of the process that persists for long 
measurement times. In the following we discuss in detail 
the behaviour of passive subdiffusive motion in terms of 
time and ensemble averages and address the peculiarities, 
that may arise for non-ergodic systems. 

Non-ergodic behaviour of the above sense is indeed ob- 
served experimentally. Fig. [2] shows the time averaged 
mean squared displacement for lipid granules in a liv- 
ing fission yeast cell. The motion is recorded by indir- 
ect tracking in an optical tweezers setup. Initially the 
granule is located in the bottom of the laser trap poten- 
tial such that the granule moves freely. Eventually the 
granule ventures away from the centre of the trap and 
experiences the Hookean trap force. As demonstrated 
in a detailed analysis the granule motion indeed exhibits 
weak ergodicity breaking, giving rise to the characteristic 
turnover from an initially linear scaling S 2 ~ A with the 
lag time A, to the power-law regime S 2 ~ A 1_a [23[. 
Moreover a pronounced trajectory-to-trajectory scatter 
is observed, again typical for systems with weak ergodi- 
city breaking. 

Free diffusion is typically quantified in terms of the 
second moment. The mean squared displacement 



J r 2 P{r,t)d 3 r 



(1) 



is obtained as the spatial average over the probability 
density function P(r, t) to find the particle at position 
r at time t. The quantity ([1]) therefore corresponds to 
the ensemble averaged second moment of the particle 
position, denoted by angular brackets, (•). In particu- 
lar, the time t enters into Eq. ([1]) only as a parameter. 
Conversely, single particle trajectories r(t) are usually 
evaluated in terms of the time averaged mean squared 
displacement defined as 



S 2 (A,T) 



1 



T- A 



T-A 



r(t + A) -r(i)) dt, (2) 



where we use an ovcrlinc ~~ ^~ to symbolise the time aver- 
age. Here A is the so-called lag time constituting a time 
window swept along the time series, and T is the over- 
all measurement time. The time averaged mean squared 
displacement thus compares the particle positions along 
the trajectory as separated by the time difference A. 




1(7 10' 10' 
A (ms) 



Figure 2: Time averaged mean squared displacement of lipid 
granule motion in fission yeast cell S.pombe, measured by op- 
tical tweezers. As function of the lag time A the initially 
linear scaling S 2 ~ A turns over to the power-law regime 
S 2 ~ A 1- ", induced by the restoring force exerted on the 
granules by the laser trap [23|. Note also the characteristic 
scatter between individual trajectories. Inset: Average beha- 
viour of the shown trajectories. In both graphs the unit of 
5 2 (A,T) is Volt 2 , i.e., the direct output voltage of the quad- 
rant photodiode. The voltage is directly proportional to the 
distance from the centre of the laser trap. 



In an ergodic system the time average of a certain 
quantity obtained from sufficiently long time series is 
equal to the corresponding ensemble mean [24], HH ■ For 
instance, for the mean squared displacement ergodicity 
would imply 



Km S 2 (A = t,T) = (r' 

T— too 



(*)>• 



(3) 



Brownian motion is ergodic, as well as certain station- 
ary processes leading to anomalous diffusion, such as 
fractional Brownian motion considered below. There ex- 
ist, however, non-ergodic processes, which are intimately 
connected to ageing properties. In what follows we dis- 
cuss two prominent models for anomalous diffusion, the 
stationary fractional Brownian motion and the ageing 
continuous time random walk, and analyse in detail the 
features of the associated time averaged mean squared 
displacement. 

Generally, anomalous diffusion denotes deviations from 
the classical linear dependence of the mean squared dis- 
placement, (r 2 (t)) ~ t. Such anomalies include ultraslow 
diffusion of the form (r 2 (t)) ~ log^ t (2(|. In contrast, 
anomalous diffusion processes may become faster than 
ballistic, for instance for systems with correlated jump 
lengths or in systems governed by generalised Langevin 
equations [27l . [28| . Here we are interested in anomalous 
diffusion with power- law dependence on time [29| , 



(v 2 (t))~2d 



r(i + a) 



(4) 
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for which the anomalous diffusion exponent belongs to 
the subdiffusive range < a < 1 , such that the limit a — 
1 corresponds to Brownian motion. The proportionality 
factor K a in Eq. (Q} is the anomalous diffusion coefficient 
of physical dimension cm 2 /sec Q . The embedding spatial 
dimension is d, and Eq. Q includes the complete Gamma 
function T(z). 

Subdiffusion of the form (QJ with < a < 1 occurs in 
the following biologically relevant systems. Fluorescently 
labelle d mRNA in E.coli bacteria cells was observed to 
follow 5 2 {A,T) ~ A a with a « 0.7 @]. This result is 
consistent with more recent findings according to which 
free RNA tracers in living cells exhibit a ~ 0.8, while 
DNA loci show a — 0.4 [ig |. Telomeres in the nucleus of 
mammalian cells were reported to follow anomalous dif- 
fusion with a 0.3 at shorter times and value a fts 0.5 at 
intermediate times, before a turnover to normal diffusion 
occurs Also larger tracer particles show anomal- 

ous diffusion, such as adeno-associated viruses of radius 
w 15 nm in a cell with a — 0.5... 0.9 |12j and endo- 
genous lipid granules, of typical size of few hundred nm 
with a « 0.75 . . .0.85 [lMl. It should be noted that 
this subdiffusion observed by single particle tracking mi- 
croscopy is consistent with results from other techniques, 
such as fluorescence correlation spectroscopy (l7H20| or 
dynamic light scattering [2l|. Subdiffusion of biopoly- 
mers larger than some 10 kD in living cells is due to 
molecular crowding, the excluded volume effect in the 
superdense cellular environment [3014321 ] . Larger tracer 
particles also experience subdiffusion due to interaction 
with the semiflexible cytoskeleton (22[. 

Knowledge of the time or ensemble averaged mean 
squared displacement of an anomalous diffusion process is 
insufficient to fully characterise the underlying stochastic 
mechanism, as the associated probability density P(r, t) 
is no longer necessarily Gaussian, and therefore no longer 
specified by the first and second moments, only [29j . This 
property is in contrast to the universal Gaussian nature 
of Brownian motion which is effected by the central limit 
theorem. At the same time the very nature of the anom- 
alous diffusion process may result in decisively different 
behaviours for diffusional mixing and diffusion-limited re- 
actions [33] . In biological cells this would imply signific- 
ant differences for signalling and regulatory processes. 
For a better understanding of the dynamics in biological 
cells and other complex fluids knowledge of the underly- 
ing stochastic mechanism is therefore imperative. 

Here we discuss the properties of continuous time 
random walk (CTRW) processes with diverging char- 
acteristic waiting time with respect to their time av- 
eraged behaviour, expanding on our earlier work [34]- 
l36j . We show that for free motion the lag time depend- 
ence of the time averaged mean squared displacement, 
(5 2 (A,T) ~ A, is insensitive to the anomalous diffusion 
exponent a for CTRW processes. In contrast, for con- 
fined CT RW subd iffusion a universal scaling behaviour 
emerges, <5 2 (A, T) ~ A 1 ^, with dynamic exponent I — a 
[curves (1") and (2") in Fig. [Q. Subdiffusion governed 



by fractional Brownian motion (FBM) leads to the scal- 
ing (5 2 (A,T) ~ A Q of the time averaged mean squared 
displacement, turning over to a saturation plateau un- 
der confinement, <5 2 (A,T) ~ A [curves (1') and (2') in 
Fig. [J. 

We particularly emphasise the irreproducible nature of 
the time averaged quantities and their associated scat- 
ter around the ensemble mean for CTRW subdiffusion 
processes. This randomness of the time averages is cap- 
tured by the distribution function of the amplitude of 
the time average. We show that even for relatively short 
trajectories this distribution is a good characteristic for 
the underlying process. In contrast for FBM processes 
the scatter typical of many single particle experiments is 
not found in the long time limit. 

In the remainder of this work, for simplicity we restrict 
the discussion to the one-dimensional case (d = I). Gen- 
eralisation to higher dimensions is straightforward. The 
article is structured as follows. We start with a brief in- 
troduction to CTRW and FBM. We then consider the 
cases of unbounded motion and confined anomalous dif- 
fusion in the subsequent two Sections. The distribution 
of the time averages will be presented thereafter. Finally, 
we discuss the velocity autocorrelation functions for sub- 
diffusive CTRW and FBM processes, before presenting a 
concluding discussion. 



II. ANOMALOUS DIFFUSION PROCESSES 

Although both CTRW and FBM give rise to an en- 
semble averaged mean squared displacement of the form 
(HJ), they are fundamentally different processes, as out- 
lined here. We note in passing that also in the random 
motion on a fractal support subdiffusion arises [33, [38| . 
We will not pursue this type of anomalous diffusion in 
the following. 



A. Continuous time random walk 

CTRW theory dates back to the work of MontroU and 
Weiss , and was championed in the analysis of charge 
carrier motion in amorphous semiconductors by Scher 
and MontroU 40]. CTRW has become a standard stat- 
istical tool to describe processes ranging from particle 
motion in actin networks [22| to the tracer motion in 
groundwater [4]| . 

Each jump of a CTRW process is characterised by a 
random jump length and a random waiting time elapsing 
before the subsequent jump. At each jump the jump 
length and waiting time are chosen independently. For 
a subdiffusive process we assume that the variance of 
the jump lengths is given by a finite value (Sx 2 ), and 
we consider the unbiased case (Sx) — 0. On a lattice of 
spacing a, we would have (Sx 2 ) = a 2 . In contrast, the 
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Figure 3: Time series x(t) for Brownian motion (top), CTRW 
(middle), and FBM (bottom). The anomalous diffusion expo- 
nent for FBM and CTRW is a = 0.5. For CTRW the stalling 
events are outstanding, while in the case of FBM strong anti- 
persistence occurs. Note that the stalling events and the anti- 
persistence are less pronounced for larger values of a. Typical 
experimental data include additional noise, such that the ap- 
pearance of measured data will not display the ideal behaviour 
shown here. 



waiting times r are drawn from the probability density 



\T(-a)\ 



(5) 



for large r, with < a < 1. This form of if> is scale-free, 
that is, the average waiting time (r) diverges, causing 
effects such as ageing [IH [43| and weak ergodicity break- 
ing [44J. In Eq. (JSJ) , the quantity r is a scaling factor. 
The anomalous diffusion constant in this case becomes 

MM 



2r Q ' 



(6) 



The scale- freeness of ip(r) allows individual waiting 
times r to become quite large. No matter how long a 
given time series is chosen, single r values may become 
of the order of the length of the entire time span covered 
in the trajectory. An example is shown in Fig. [3j the 
stalling events with large waiting times r are quite dis- 
tinct. For values of the anomalous diffusion exponent 
a that are closer to 1 the stalling is less pronounced. 
Physically, the power-law form of the waiting time dis- 
tribution ip(r) may be related to comb models 37[ or 
random energy landscapes [H,|46j]; for more details, com- 
pare Refs. 0413. 



B. Fractional Brownian motion 

FBM is a Gaussian process with stationary increments. 
Its position is defined in terms of the Langevin equation 



dx(t) 
dt 



or, alternatively, 



x(t) = / i{t')dt' 
Jo 



(J) 



(8) 



The motion is driven by stationary, fractional Gaussian 
noise with zero mean (£(t)) and long-ranged noise 
correlation [FlJ 

<£(fi)£(ta)) = aK*(a-l)\h-t2\ a - 2 

+2aK*\t 1 -t 2 \ a - 1 6{t 1 -t 2 ), (9) 

contrasting the uncorrelated noise for normal diffusion 
a = 1: (£(ii)£(i 2 )} = WiSfa - t 2 ). In Eq. © the 
anomalous diffusion exponent is connected to the tra- 
ditionally used Hurst exponent by H = a/2, and we 
introduce the abbreviation K* a = K a /T(l + a) for con- 
sistency with standard FBM notation. For subdiffusion 
the fractional Gaussian noise is anticorrelated, decaying 

like (f(ti)£(i 2 )) K * a a\a - l\\ti - i 2 | Q " 2 - This implies 

that a given step is likely to go into the direction oppos- 
ite to the previous step. The corresponding oscillating 
behaviour is seen in Fig. |3l The position autocorrelation 
function of FBM becomes 



(x(ti)x(t 2 ))=K*Jt? + t%-\t 1 -t 



(*? 



(10) 



so that at equal times t\ = t 2 we recover the mean 
squared displacement Q. 

If the fractional Gaussian noise is not considered ex- 
ternal, but the validity of the fluctuation-dissipation the- 
orem is imposed, one obtains the generalised Langevin 
equation (GLE) |52| . 



d 2 x(t) 
1 dt 2 



+ vt(t), (11) 



where we write /3 instead of the exponent a in Eq. ©. 
In the GLE (|III) . we defined the coupling coefficient 

r\ = ^JkBTj/[(3Kp((3 — 1)] according to the fluctuation- 
dissipation theorem, where kg is the Boltzmann constant 
and T the absolute temperature. In what follows we only 
consider the overdamped limit, in which the inertia term 
md 2 x(t)/dt 2 can be neglected. The GLE then gives rise 
to the form 



2-/3 



(12) 



of the mean squared displacement. In contrast to FBM, 
that is, the GLE leads to subdiffusion for persistent noise 
with 1 < f3 < 2, while < (3 < 1 yields superdiffusion. 
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FBM and the related GLE are used to describe pro- 
cesses such as long term storage cap acity of water reser- 
voirs climate fluctuations (54 |, e conomical market 
dynamics (55[, single file diffusion 56|, and elastic mod- 
els [I?) ■ Motion of this type has also been associated 
with the relative motion of aminoacids in proteins (58l |. 
and the free diffusion of biopolymers under molecular 
crowding conditions [lfj, lljj [59j . 



III. FREE ANOMALOUS DIFFUSION 

Let us begin with considering anomalous diffusion on 
an infinite domain and without drift. To find an ana- 
lytical expression for the time averaged mean squared 
displacement @ we note that even a Brownian process 
recorded over a finite time span T will show fluctuations 
in the number of jumps performed during T. To average 
out these trajectory-to-trajectory fluctuations we intro- 
duce the ensemble mean, 



S 2 (A,T) = 



1 



T — A 



T-A 



(x{t + A) - x(t) 



dt. 



We can then express the integrand in terms of the vari- 
ance of the jump lengths, (Sx 2 ), and the average number 
of jumps n(t,t + A) in the time interval {t,t + A), as 
follows: 

({x{t + A) -z(<)) 2 ^ = (5x 2 )n(t,t + A). (14) 

For a regular random walk on average every jump occurs 
after the waiting time (r). Thus n(t,t + A) — A/(r), 
and 



8^{A,T)) = 2K l A, 



(15) 



where we defined the diffusion constant K\ = 
(5x 2 ) /[2(t)]. In the Brownian limit the time averaged 
mean squared displacement (|15p in terms of the lag time 
A takes on exactly the same form as the ensemble aver- 
aged mean squared displacement ([T]) as function of time t. 
This is not surprising, as Brownian motion is an ergodic 
process. For sufficient duration T of the time records any 
time average converges to the corresponding ensemble av- 
erage, and thus the ensemble average in expression (fT5|) 
is no longer necessary. 



A. Continuous time random walk 

The number of jumps of a CTRW process with waiting 
time distribution of the form ([5]) on average grows sublin- 
early with time, n(0,t) ~ t a / [r a T (1 + a)} '47]. This time 
evolution translates into the mean squared displacement 
((T|) with the anomalous diffusion coefficient ^ . Combin- 
ing the time dependence of n(0, t) for CTRW subdiffusion 




Figure 4: Time averaged mean squared displacement for un- 
confmed CTRW motion with a — 0.5, shown for 20 individual 
trajectories. The overall measurement time is T = 10 (a.u.). 
Note the local changes of the slope in the trajectories as well 
as the complete stalling in the lowest curve, all bearing wit- 
ness to the scale-free nature of the underlying waiting time 



(13) distribution i[>(t) ~ t q /t 



with the definition (fT3"| we obtain the following result for 
the time averaged mean squared displacement, 



<5 2 (A,T)) ~2K a 



(16) 



in the limit A <C T [3J, |6fJ • This follows from expansion 
of the relation n(t,t + A) = n(0,t + A) - n(0,t). The 
noteworthy feature of this result is that the linear lag 
time dependence of normal diffusion remains completely 
unaffected by the anomalous nature of the stochastic pro- 
cess. Only the dependence on the overall measurement 
time T witnesses the underlying subdiffusion. Eq. (|16[) 
can be understood if we notice that for normal diffusion 
we have S 2 (A,T) ~ A/(t) [since according to Einstein 
K\ oc 1/(t}] where (r) is the mean time between jump 
events. Now, for the subdiffusive case (t) diverges and 
must be replaced by iP(t)tcIt oc T 1 ^ , which explains 
the term S 2 oc A/T 1_Q . Fig. [4] shows 20 trajectories of 
a CT RW proc ess with a — 0.5, displaying the general 
trend S 2 (A,T) ~ A. 

Rewriting the result (|16l) in the form 



(17) 



<5 2 (A,T)) - 2K(T)A 



we see that the effective diffusion coefficient K(T) de- 
cays as function of the measurement time T. The longer 
the system evolves after its initial preparation, the less 
mobile it appears, consistent with the ageing property of 
CTRW subdiffusion. For instance, in the picture of the 
trapping events, in the course of time deeper and deeper 
traps may be encountered by the diffusing particle, such 
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that it gets more and more stuck. This increasing immob- 
ility is mirrored in the form K(T). Note that Eqs. (fT6|) 
and (JXTJ) suggest that from measuring the A dependence 
of an anomalous diffusion process one might draw the er- 
roneous conclusion that normal diffusion were observed. 

The disparity between the time averaged mean squared 
displacement (1161) and its ensemble averaged counterpart 
(|4j) demonstrates the weak ergodicity breaking charac- 
teristic of a process with diverging time scale [IH, 13]. 
Single or few, long waiting times are also responsible 
for the pronounced deviations between individual realisa- 
tions shown in Fig. 21 This apparent irreproducibility of 
the time averaged mean squared displacement is again in- 
timately coupled to the weak ergodicity breaking nature 
of the CTRW subdiffusion process. We will discuss this 
feature more quantitatively in terms of the distribution 

4> a (0 as function of the relative deviation £ = S 2 / ^<5 2 ^} 

of the time averaged mean squared displacement around 
its ensemble mean in Section [V] 
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Figure 5: Time averaged mean squared displacement for un- 
confined FBM with a = 0.5. The scatter between individual 
trajectories is minute, mirroring the ergodicity of this process. 



B. Fractional Brownian motion 

In contrast to the above behaviour of CTRW subdif- 
fusion processes, FBM is ergodic. Indeed, by help of 
the position autocorrelation (|10[) the time averaged mean 
squared displacement (TT51) becomes 61 1 



5 2 (A,T)) = 2iY*A c 



(18) 



exactly matching the ensemble averaged mean squared 
displacement, (x 2 (t)) = 2K*t a . However, ergodicity is 



reached algebraically slowly [61], see also below. In Fig. [5] 
the comparatively minute scatter between individual tra- 
jectories supports the ergodic behaviour of FBM pro- 
cesses. The somewhat larger deviations at longer lag 
times are due to worsening statistics when A — » T. 



IV. CONFINED ANOMALOUS DIFFUSION 

Since the motion of tracer particles is typically con- 
fined, for example, by the cell walls or internal mem- 
branes, we now consider the important case of anomalous 
diffusion in a bounded domain. 

In a finite interval [-L, L] the mean squared displace- 
ment of a Brownian particle initially released well away 
from the boundaries will grow linearly in time and even- 
tually turn over to a stationary plateau of magnitude 
(x 2 ) s t — L 2 /3. Similarly, if the particle evolves in the 
confinement of an external harmonic potential of the form 
V(x) = ^muj 2 x 2 , the thermal value (ir 2 }th = fcsT/[mw 2 ] 
will eventually be attained. 



A. Continuous time random walk 

How is this behaviour modified for a CTRW subdif- 
fusion process? Under the influence of an arbitrary ex- 
ternal potential V(x) = — f x F(x')dx' defining the force 
F(x), CTRW subdiffusion can be described by the frac- 
tional Fokker-Planck equation [29|, or, equivalently, 
in terms of the following coupled Langevin equations (62j | : 



dx(s) 

ds 
dt(s) 

ds 



K 
u){s) 



F(x) + r)(s) 



(19a) 
(19b) 



Here the position x is expressed in terms of the parameter 
s (the internal time), and driven by the white Gaussian 
noise rj(s). Thus, Eq. (|19al) defines standard Brownian 
motion x(s), where K is the diffusivity for the normal 
diffusion process in internal time s. Laboratory time t 
is introduced by the so-called subordination through the 
process w(s), given by the probability density function 



Pt{s) 



K 



l/c 



t 



,•1+1/" 



( 


~K a - 


1/a t \ 




K 


s x+i/ a J 



(20) 

where l a (z) is a one-sided Levy stable probability dens- 
ity with Laplace transform J Q l a {z) exp(—uz)dz = 
exp(— u a ). Thus, Eq. (I19b[) transforms the Brownian 
process x(s) with diffusivity K into the subdiffusive mo- 
tion with generalised diffusivity K a . On the basis of this 
scheme, the ensemble averaged position-position correla- 
tion in an arbitrary confining potential V(x) becomes 



{X(ti)x{t2)) = (V)i 



B(ti/t 2 ,a, 1 - a) 
r(a)r(l - a) 



(21) 
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valid in the limit t 2 > h > (1/^Ai]) 1 /", where Ai 
is the smallest non-zero eig envalue of the corresponding 
Fokker-Planck operator [35]]. Eq. (|2T|) demonstrates that, 
despite the confinement, the process is non-stationary. In 
result (|2T|) we used the incomplete Beta function 



B(ti/t 2 , a, 1 - a) = 



*l/*2 



z a_1 (l - z)- a dz (22) 



and defined the first and second Boltzmann moments, 
whose general definition is 



f f°° 

(x n )B = ^ y_ 00 a:nexp 

The partition function reads 



k B T 



dx. 



2? = 



exp 



k B T 



dx. 



(23) 



(24) 



At large time separation, t 2 3> ti, the position autocor- 
relation decays algebraically, 



(x(h)x(t 2 )) 



aT(l + a)T(l - a) 



+ <*>!, (25) 



towards the value (x) B . The corresponding limiting be- 
haviour of the incomplete Beta function reads 



B(t/(t + A),a,l-a) 
r(a)r(l-a) 



1 - 



sin(7ra!) / A 



(1 - a)Tt \t 



l-a 



(26) 



Inserting expression (|2"Tj) into the definition of the time 
averaged mean squared displacement ([TB")) we obtain [3o| 



S 2 (A,T) 



T 



T-A 



(x(i + A)x(t + A)) 



+(x(t)x(i)) - 2(x(t + A)x(i)) 



rff. 



(27) 



Then, with relation (f26|) we arrive at the scaling beha- 
viour 



^(A,T))~((x 2 ) B -(x>|) 



2 sin(7ra) 
(1 — a)ira 



(28) 

valid in the limits A/T « 1 and A » (l/fA^Ai]) 1 /". 
Result (|28p is quite remarkable: instead of the naively as- 
sumed saturation toward the stationary plateau a power- 

(A/T) 1 -" is observed. Only 



law growth <S 2 (A,T) 

when the lag time A approaches the overall measurement 
time T the singularity in expression (|13j) causes a dip to- 
ward the plateau of the ensemble averaged mean squared 
displacement. Additionally Eqs. (|2"Tj) and (|2"8j) are uni- 
versal in the sense that the exact form of the confining 
potential solely enters into the prefactor through the first 
two moments of the Boltzmann distribution correspond- 
ing to the confining potential V(x). We note that the 
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A 



Figure 6: Simulated behaviour of ^<5 2 (A)y for an harmonic 

binding potential V(x) — 2x 2 (□) and a particle in a box 
of length 2 (O); with a = 1/2, measurement time T = 10 7 
(a.u.), fc B T = 0.1, and K .s = 0.0892. Without fitting, the 
lines show the analytic results for the transition from ~ A 1 

for short lag times according to Eq. (|16[) ( and — ■ — ) 

to ~ A 1_Q for long lag times, Eq. ([28]) (— ). For long A, 

^<5 2 (A)^ /(x 2 )b exhibits universal behaviour, in the sense 

that the curve does not depend on the external field. 



asymptotic scaling ~ (A/T) 1 a is consistent with the 
numerical analysis presented in Ref. (631 ). 



Fig. [6] depicts the ensemble mean of the time averaged 
mean squared displacement (|13p for two types of confine- 
ment, an harmonic potential and a box potential. The 
particle is initially placed at the bottom of the poten- 
tial well and in the middle of the box potential, respect- 
ively. This is why at short lag times the particle exhibits 

A typical for unconfined 



the linear scaling ^<5 2 (A,T) 

motion, before turning over to the confinement-induced 
scaling ~ A 1_Q . Note that this plot is fit-free, i.e., the 
theoretically calculated asymptotic behaviours nicely fall 
on top of the simulations results. 



As shown in Fig. [7] individual trajectories still ex- 
hibit the pronounced scatter typical for CTRW subdiffu- 
sion. Visually the scatter does not change between the 
unbiased initial motion and the confinement-dominated 
part of the process. This is due to the fact that the scat- 
ter is caused by the scalc-frccness of the waiting times. 
In our process waiting times and jump lengths are de- 
coupled, corresponding to the subordination property of 
CTRW subdiffusion [29|. 
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Figure 7: Scatter between individual trajectories in the 
harmonic potential V(x) = 2x 2 . The extent of the de- 
viations between the trajectories does not change qualit- 
atively between the initial unbiased motion and the later 
confinement-dominated regime. The squares (□) repres- 
ent simulations results for the ensemble average (s 2 (A,T) 
Same parameters as in Fig. 



B. Continuous time random walk with waiting 
time cutoff 



What happens if we introduce a cutoff in the power-law 
of the waiting time distribution of the form 



d 



1 



(r + r)' 



exp 



r/r*) 



(29) 



in which a characteristic time scale r* is introduced, ter- 
minating the power-law scaling? For the waiting time 
density (|29|) the characteristic waiting time / °° Tip{r)dT 
becomes finite. At sufficiently short times one would ex- 
pect this process to still exhibit the features of CTRW 
subdiffusion, while at times r 3> r* the process should 
converge to regular Brownian motion with a Gaussian 
propagator. In Fig. [5] we demonstrate that for a suitable 
choice of the cutoff time r* the beha viour of t he time av- 
eraged mean squared displacement <5 2 (A,T) under con- 
finement preserves the characteristic non-ergodic features 
of CTRW subdiffusion, i.e., the turnover from the initial 
scaling ~ A to the confinement-dominated scaling A 1- ". 
Below we will show that at the same time the distribu- 
tion of the time average around its ensemble mean is 
significantly altered. 



C. Fractional Brownian motion 

Ergodicity remains unaffected for FBM, that is con- 
fined to an interval of size [—L,L]. Namely, conver- 
gence of <5 2 (A, T) is observed toward the stationary value 



Figure 8: CTRW subdiffusion for a waiting time with ex- 
ponential cutoff, Eq. ((29} , in the box [—3, 3] with reflecting 
boundary conditions. The generic behaviour of initial and 
final scaling, <5 2 (A,T) ~ A and ~ A 1 "" remains unaltered. 
We chose a = 0.8, t = 1, r* = 20 and overall measurement 
time T = 1000 (a.u.). 



(a; 2 } st = L 2 /3 64]. For FBM under the influence of 



an harmonic potential V(x) 



,2„2 



the position- 



position correlation can be calculated exactly [65|. In- 
serting into the ensemble mean (|13p of the time averaged 
mean squared displacement one can show that the ini- 
tial behaviour (x 2 (t)) ~ t a turns over to the stationary 
value (x 2 ) st = r(a + l)fcsT/[mw 2 ]. In general, the time 
averaged mean squared displacement for confined FBM 
converges to a constant: 



S 2 (A,T)) ~ const. 



(30) 



V. FLUCTUATIONS OF THE TIME AVERAGE 
AND ERGODICITY BREAKING PARAMETER 

The devi ations of the time averaged mean squared dis- 
placement (5 2 (A,T) between individual trajectories can 
be quantified in terms of the probability density function 
4> a {C) of the dimensionless ratio 



S 2 (A,T) 



5 2 (A,T) 



(31) 



of the time averaged mean squared displacement over its 
ensemble mean. For an ergodic process this distribution 
is necessarily sharp, 



0erg(£) = - 1), 



(32) 



for long measurement time T. Deviations from this form 
are expected for non-ergodic processes such as CTRW 
subdiffusion, but also for relatively short trajectories. We 
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here address both effects and demonstrate that the distri- 
bution <j) a (£) is a quite reliable means to distinguish dif- 
ferent stochastic processes even when the recorded time 
series are fairly short. 



A. Continuous time random walk 

For CTRW subdiffusion with power-law waiting time 
distribution ([5]) the distribution assumes the form [34], HB| 



MO 



r(l + q) 1 /" / T(l + a ) 

ofl+l/a la { £l/o 



\l/a 



(33) 



for T — > oo, where / a (z) is again a one-sided Levy stable 
distribution, whose Laplace transform is J£{l a (z)} — 
exp (— u a ). Special cases include the Brownian limit (j3"2")l 
for a = 1 and the Gaussian shape 



'1/2 



(0 



■ exp 



e 



(34) 



for a = 1/2. In fact, we observe an exponentially fast 
decay of </> a (£) with large £ for all < a < 1, compare 
Appendix [Al 

For short trajectories the basic shape of the distribu- 
tion (f>a((,) is surprisingly well preserved [67^ . This is 
demonstrated in Fig. El in which the characteristic asym- 
metric shape with respect to the ergodic value £ = 1 for 
the case a = 0.5 is reproduced for a process with T = 128, 
even for a lag time as long as A = 100. No significant de- 
pendence on the confinement is observed. Also for larger 
values of a the finite value for small £ is similarly re- 
produced for short trajectories [63] , pointing at the quite 
remarkably reliability on the shape of the scatter distri- 
bution <fi a (£) for CTRW subdiffusion. As we show below, 
the distribution for FBM processes with its zero value at 
£ = can be clearly distinguished from the CTRW form. 

A useful parameter to quantify the violation of ergodi- 
city is the ergodicity breaking parameter [3~i| 



EB = lim 



(F) ; 



-(s 2 



2r(i 



— \ 2 

s 2 



r(l + 2a) 



1 (35) 



EB varies from EB = 1 for a — > monotonically to 
EB = in the Brownian limit a = 1. For the special 
case a — 1/2 one finds EB = ir/2 — 1 rs 0.57, while for 
a — 0.75, EB « 0.27. 



B. Continuous time random walk with waiting 
time cutoff 



The scatter distribution is sensitive to few extreme 
events. When these are lacking, the distribution should 
be significantly different from the form discussed for 
CTRW subdiffusion with diverging characteristic time 
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Figure 9: Distribution 4>a(£,) of the time averaged mean 
squared displacement, with £ = S 2 / (s 2 ^ for a CTRW pro- 
cess with a — 0.5, r = 1, and T = 128 (a.u.). The filled and 
open symbols, respectively, represent the unconfined case and 
confined motion on the interval [—2,2]. We see quite good 
agreement with the expected Gaussian limiting distribution 
of the scatter, Eq. (|34p . centred around £ = 0. In the inset 
we show the case for the largest measured lag time, A = 100, 
also in good agreement with the predicted shape of the dis- 
tribution. 



scale. Indeed, as shown in Fig. [TU] for the cutoff waiting 
time distribution defined through Eq. (f29)) . drops 
down to zero around £ = and assumes an almost Gaus- 
sian shape around the ergodic value £ = 1. Note that the 
simulations were performed with the same parameters as 
for Fig. [51 in which the breaking of ergodicity still per- 
sists, despite the cutoff. This demonstrates that different 
quantities have different sensitivity to the cutoff. 



C. Fractional Brownian motion 

FBM is based on long-ranged correlations. However, 
we can obtain an approximate expression for the scatter 
distribution K37fl 



T 



47TT-t 



■ exp 



(£-i) 2 (r-A) 

4 r t 



(36) 



where is an intrinsic time scale. Expression (|36[) is 
valid for sufficiently small lag times A, for which the 
correlations are neglected [67j ■ Note that this distribu- 
tion is independent of a and centred around the ergodic 
value £ = 1. In particular, in the long measurement time 
limit T — > oo, the distribution converges to the sharp 
form 4>(£) = <5(£ — 1). Also this behaviour is surprisingly 
well preserved for short trajectories, as demonstrated in 
Fig. HU 
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Figure 10: Scatter distribution 4> a (0 for CTRW subdiffusion 
with cutoff, as defined in Eq. (|29p . The observed form is 
clearly different from the much more asymmetric shape in 
Fig. [9] The value of <j> a (£) reaches zero for small values of £• 
An almost Gaussian shape centred around the ergodic value 
£ = 1 is observed. Same parameters as in Fig. [8] 
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Figure 11: Scatter distribution of the time averaged 

mean squared displacement, with £ = S 2 / (^5 2 ^ for FBM sub- 
diffusion with a = 0.5 and T — 128 (a.u.). The filled and open 
symbols, respectively, represent the cases for unconfined and 
confined motion. For small values of the lag time A we see 
quite good agreement with the expected Gaussian limiting 
distribution of the scatter, Eq. (|36p . that is centred around 
£ = 1. For larger values of A deviations are observed, how- 
ever, even for A = 100 the value around £ = is consiste ntly 
zero. A similar behaviour is found for larger values of a [67| . 
Note that when T — > oo the process is ergodic and <fi(£) ap- 
proaches a 8 function. 



For FBM one can define a quantity similar to the er- 
godicity breaking parameter (|55|) . Namely, without tak- 
ing the long measurement time limit, we obtain the nor- 
malised variance of the time averaged mean squared dis- 
placement 



V 



<5 2 (A,T) 



S 2 (A,T) 



<$ 2 (A,T) 



(37) 



It turns out that the convergence to ergodicity is algeb- 
raically slow, and for subdiffusion < a < 1 one obtains 
a decay, which is inversely proportional to T 61]: 



The coefficient k(a) here is defined as [61 1 

k(a) = J ((i + l) Q + |£-l| Q -2t Q )di. 



(38) 



(39) 



It increases continuously from zero at a = to k(l) = 
4/3. 

The quantity V is connected to the probability density 
(|36) by V = 2t7(T - A) - obtained from ap- 

proximating that the values of £ may range in the interval 
(—oo, oo), which appears reasonable given the sharp de- 
cay of 4>(£) to at £ = for FBM. Apart from the fact 
that V explicitly depends on a while the approximation 
leading to Eq. (pio 7 )) loses the a dependence, both quant- 
ities decay ~ 1/T, and the internal time scale takes 
on a role similar to the lag time A. 



VI. VELOCITY AUTOCORRELATION 
FUNCTIONS 

A typical quantity accessible from experimental data 
is the velocity autocorrelation function, which is defined 
through 



Cl e) (r) = ^( Kt + e) - x(t)) (x(e) - x(0)) ). 



(40) 



Here the velocity is defined as the difference quotient 
v(t) = e _1 [a;(r + e) —x(t)]. The velocity autocorrelation 
(P4U1) was suggested as a tool to distinguish between dif- 
ferent subdiffusion models [10(. We show here that for 
confined processes the shape of the velocity autocorrela- 
tion function does not allow for a significant distinction 
between subdiffusive CTRW and FBM. Note that the cal- 
culation of Ci e \r) amounts to determine four two-point 
correlation functions of the type (x(ti)x(t2)) ■ 



A. Continuous time random walk 

For unbounded CTRW process starting with initial 
condition a;(0) = 0, the position correlation function is 
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Figure 12: Normalised velocity autocorrelation function 
Ci e) (T) for an unconfined CTRW process with a = 0.5. The 
simulations (O) were performed for time 10000 (a.u) with 
e — 10 (a.u). The theoretical behaviour ( — ) is given by 
Eq. g2j). 



(41) 



given by the following expression [68( , 

2K n 



(x(h)x(t 2 )) 



r(i + a) 



min{ti,i 2 } 



This result for free CTRW is due to the fact that the 
jump lengths in the interval (£1,^2) for t 2 > t\ have zero 
mean. With Eq. (|41[) we find 



d 6) (r) 







t < e 
t > € 



(42) 



The velocity autocorrelation function for an unbounded 
CTRW process does not yield negative values due to the 
absence of correlations for different jumps. In this case 
the velocity autocorrelations can be easily distinguished 
from the behaviour of FBM processes, see below. Note 
that Cv (t) in Eq. (|4"2"j) is non-analytic at r = e, an ob- 
servation that is confirmed in simulations. The behaviour 
of the velocity autocorrelations for CTRW subdiffusion 
are displayed in Fig. Q2] 

For a confined CTRW process the situation is quite dif- 
ferent. To explore the behaviour of the corresponding ve- 
locity autocorrelation function we use the general result 
for the correlation function of confined CTRW, Eq. (|2T| . 
For simplicity we assume a symmetric potential such that 
(x)b = 0. With the initial condition x(0) =0we obtain 



1 



CW( ) T(a)T(l-a) 
x<B ^l-a)- fl ( W - )e>r 

Figs. [T3] and Q3] (for the absolute value) display excel- 
lent agreement of Eq. (|4"5|) with simulations of CTRW 



Figure 13: Normalised velocity autocorrelation function 
C^(r) for a free FBM process and a confined CTRW process 
with a — 0.5. The CTRW simulations (O) were performed 
over the time range 10000 (a.u), and the system corresponds 
to a particle on a lattice with ten lattice points and reflecting 
boundaries. We chose e = 10 (a.u). The theoretical predic- 
tion for the CTRW (— ) is given by Eq. fl43j), and for FBM (- 
-) by Eq. ([45]). 



subdiffusion for a lattice of size 10 with reflecting bound- 
ary conditions, and a — 1/2. We observe that for con- 
fined motion the CTRW velocity autocorrelation func- 
tion indeed attains negative values and has a minimum 
on t = e, as the confinement effectively induces correl- 
ations. For long t the velocity autocorrelation function 
decays to zero (from the negative side) as the power-law 



d e) (o) 

valid for t > e. 



1 



r(a)r(i - a) \t 



(r) ! 



(44) 



B. Fractional Brownian Motion 

For free FBM the position correlation function (flU)) 
together with the definition ([40]) produce the result 



d £) (r) 

ci e) (o) 



(t + e)° - 2t" +jr-ej 
2e a 



(45) 



This function yields negative values for sufficiently long 
t, and its minimum value (2" _1 — 1) is assumed at r = e. 
For long r it decays toward zero from the negative side 
in the power-law form 



c4 £) (t) 



2-a 



(46) 



valid for r 3> e. 

We see that both the velocity autocorrelation function 
of the confined CTRW and the one for free FBM acquire 
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Figure 14: Absolute value of the normalised velocity auto- 
correlation function (r) for confined CTRW process with 
a = 0.5. The CTRW simulations (O) are based on the para- 
meters from Fig. 1131 The theory line ( — ) corresponds to 
Eq. (|43l). 
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Figure 15: Normalised velocity autocorrelation function 
Ci e) (r) for confined FBM with a = 0.5. The simulations 
were performed on the interval [—2,2], and e = 10 (a.u.). 
For comparison, the velocity autocorrelation for free FBM is 
drawn based on Eq. (I45p . 



negative values and decay as power-laws for large r, see 
below. Further both obtain a sharp minimum for r = 
e. Confined FBM behaves similarly to free FBM, see 
Fig. [15j From our analysis it becomes clear that the 
shape of the velocity autocorrelation function may not be 
a good diagnosis tool to distinguish between subdiffusivc 
CTRW and FBM processes. We note that for CTRW 
subdiffusion the time averaged correlation functions will 
be random variables, exactly like the time averaged mean 
squared displacement. 



VII. DISCUSSION 

The physical mechanisms leading to subdiffusion in 
biological cells or other complex systems are varied. 
The anomalous diffusion of inert biopolymers larger than 
some lOkD in living cells is due to molecular crowding, 
these are excluded volume effects in the viscoelastic, su- 
perdense cellular environment (30M32^ . Potentially sim- 
ilar effects occur in single file diffusion: A tracer particle 
diffusing in a one dimensional channel and interacting 
with other Brownian particles will slow down its ran- 
dom motion, and subdiffusion with a = 1/2 is observed 
What is the resulting motion like? Consider a 
particle, that gets repeatedly trapped during its random 
walk. Such trapping may give rise to a broad distribu- 
tion of waiting times, as embodied in the CTRW model. 
Such broad distribution of trapping times could be due 
to chemical binding of the tracer particle to its environ- 
ment on varied time scales, or due to active gelation/de- 
gelation of the crowding particles. Long-tailed trapping 
time distributions were observed for /mi-sized particles 
due to interaction with a semiflexible actin mesh (22j. 
There, the anomalous diffusion exponent a depends on 
the size of the tracer particles versus the typical actin 
mesh size. Alternatively, motion patterns of FBM type 
may be due to coupling to the viscoelastic crowding envir- 
onment [l!| . Yet other approaches are based on polymer 
dynamics. De Gennes' reptation model of a polymer dif- 
fusing in a tube formed by foreign chains in a melt yields 
subdiffusion with a = 1 /4, and was used as a possible ex- 
planation of the short time dynamics of telomere motion 

One may speculate whether the widely observed sub- 
diffusion of biomolecules and passive tracers is a mere 
coincidence of nature that should be attributed to the 
dense environment in the cell and/or to the wide distri- 
bution of obstacles and reactions in the cell. Or, maybe, 
subdiffusion is by itself a goal which was obtained via 
evolution? There exist claims that sub-diffusion is use- 
ful in certain cellular search strategies [j| |(59| . One may 
also view subdiffusion to be a sufficiently slow process to 
help maintaining the organisation of the genome in the 
nucleus of the cell without the need for physical com- 
partments pr| . Since reactions in many cases are con- 
trolled by diffusion, the emergence of slower-than-normal 
diffusion has far-reaching implications to signalling and 
regulatory processes in the cell. The standard theor- 
ies of diffusion-controlled reactions under such circum- 
stances must be replaced by subdiffusion-controlled re- 
action models. While these issues are clearly important, 
we here focused more on the characterisation of traject- 
ories of single particles in the cell, in particular, as a 
quantitative diagnosis method to probe the nature of the 
stochastic motion in cells and other complex media. 

Single particle tracking microscopy is a widely used 
technique. It allows one to locally probe complex systems 
in the liquid phase in situ. In particular it has become 
one of the standard tools in biophysical, colloidal, poly- 
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meric, and gel-like environments. The motion of tracers 
in these systems is often anomalous. Typically the exper- 
imentally recorded time series are evaluated in terms o f 
the time averaged mean squared displacement <5 2 (A,T). 
Here we collected the behaviour of 8 2 (A,T) for the two 
most prominent anomalous stochastic processes, the con- 
tinuous time random walk and fractional Brownian mo- 
tion. 

For CTRW subdiffusion, connected with ageing and 
weak ergodicity breaking, in an unconfined system the 
time ave raged me an squared displacement scales linearly 
in time, <5 2 (A,T) [Eq. (|16p ]. renouncing the anomalous 
nature of the process. Only the dependence on the over- 
all measurement time T pays tribute to the underlying 
subdifFusion. This behaviour contrasts the anomalous 
scaling (x 2 (t)) = 2K a t a of the ensemble averaged mean 
squared displacement. Under confinement no plateau 
value is observed as in the ensemble average, instead, 
a power-law of the form S 2 (A,T) ~ A 1 "", Eq. flM}, is 
found. Interestingly, this characteristic behaviour is ap- 
proximately preserved when an appropriate cutoff in the 
waiting time is introduced. Conversely the distribution 
of the time averaged mean squared displacement around 
its ensemble mean becomes almost Gaussian in presence 
of the cutoff, while it has an exponential decay and a fi- 
nite value at £ = when the system is non-ergodic and 
exhibits ageing. 

FBM, in contrast, is ergodic, although ergodicity is 
reache d algebr aically slowly. Free FBM subdiffusion 
shows (5 2 (A,T) ~ A Q , Eq. (fLS)) . while under confinement 
the plateau value of the ensemble average is reached, 
Eq. (|3"0)) . For finite trajectories the distribution of the 
time averaged mean squared displacement is approxim- 
ately Gaussian for short lag times. 

Our analysis demonstrates how different the two 
stochastic processes are, despite sharing the same form 
of the ensemble averaged mean squared displacement. 
At the same time the velocity autocorrelation of con- 
fined CTRW subdiffusion is hardly distinguishable from 
that of FBM subdiffusion. Given a recorded time series 



of anomalous diffusion from experiment it is important 
to know more precisely which stochastic process is re- 
sponsible for the observed behaviour, in particular, with 
respect to diffusion-limited reactions, general transport 
behaviour, and related processes such as gene regulation. 
The analysis presented here, al ong with complementary 
tools discussed in Refs. [23, [59|, [7Q| , will be instrumental 
in the classification of anomalous diffusion behaviour. 
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Appendix A: Scatter distribution 

Using the theory of Fox //-functions we obtain the ex- 
act form for the distribution </> Q (£) of the dimensionless 
variable £ [zj: 



= 



l/c 



(0,1) 
(0,1/a) 



This function has the series expansion 



n=0 

and the asymptotic form 

(l-a)/(2a)-a 



n!r(— an) 



(Al) 



(A2) 




V(i-") N 



(A3) 
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